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Quantum phases of matter are characterized by the underlying correlations of the many-body 
system. Although this is typically captured by a local order parameter, it has been shown that 
a broad class of many-body systems possesses a hidden non-local order. In the case of bosonic 
Mott insulators, the ground state properties are governed by quantum fluctuations in the form of 
correlated particle-hole pairs that lead to the emergence of a non-local string order in one dimension. 
Using high-resolution imaging of low-dimensional quantum gases in an optical lattice, we directly 
detect these pairs with single-site and single-particle sensitivity and observe string order in the 
one-dimensional case. 



The realization of strongly correlated quantum many- 
body systems using ultracold atoms has enabled the di- 
rect observation and control of fundamental quantum ef- 
fects [TH3]. A prominent example is the transition from 
a superfluid (SF) to a Mott insulator (MI), occurring 
when interactions between bosonic particles on a lattice 
dominate over their kinetic energy [4-8 . At zero tem- 
perature, and in the limit where the ratio of kinetic en- 
ergy over interaction energy vanishes, particle fluctua- 
tions are completely suppressed and the lattice sites are 
occupied by an integer number of particles. However, at 
a finite tunnel coupling, but still in the Mott insulating 
regime, quantum fluctuations create correlated particle- 
hole pairs on top of this fixed-density background, which 
can be understood as virtual excitations. These particle- 
hole pairs fundamentally determine the properties of the 
Mott insulator such as its residual phase coherence [9] 
and lie at the heart of superexchange-mediated spin in- 
teractions that form the basis of quantum magnetism in 
multi-component quantum gas mixtures |10H12j . 

In a one-dimensional system, the appearance of cor- 
related particle-hole pairs at the transition point from 
a superfluid to a Mott insulator is intimately connected 
to the emergence of a hidden string-order parameter Op 




Here 5hj = fij — h denotes the deviation in occupation 
of the jth lattice site from the average background den- 
sity, and k is an arbitrary position along the chain. In 
the simplest case of a Mott insulator with unity filling 
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(fl = 1), relevant to our experiments, each factor in the 
product of operators in Eq.[T] yields —1 instead of +1, 
when a single-particle fluctuation from the unit back- 
ground density is encountered. In the superfluid, particle 
and hole fluctuations occur independently and are uncor- 
rected, such that Op = 0. However, in the Mott insu- 
lating phase, density fluctuations always occur as cor- 
related particle-hole pairs, resulting in Op ^ 0. For a 
homogeneous system, Op is expected to follow a scal- 
ing of Berezinskii-Kosterlitz-Thouless (BKT) type |15j . 
Non-local correlation functions, like the string-order pa- 
rameter defined above, have been introduced in the con- 
text of low-dimensional quantum systems. They classify 
many-body quantum phases that are not amenable to 
a description through a local order parameter, typically 
used in the Landau paradigm of phase transitions. Exam- 
ples include spin-1 chains [16] and spin- 1/2 ladders [17] . 
fermionic Mott and band insulators [TS] and Haldane in- 
sulators in one-dimensional Bose gases [131 HI]. Recently, 
the intimate connection of string order and local symme- 
tries has been uncovered [H] and wide-ranging classifi- 
cation schemes for quantum phases using such symmetry 
principles have been introduced [501 121] ■ Here we show 
for the first time that correlated particle-hole pairs and 
string order can be directly detected using single-atom- 
resolved images of strongly correlated ultracold quantum 
gases [22] 123]. 

We prepared a two-dimensional degenerate gas of ul- 
tracold 87 Rb atoms, before shining in a two-dimensional 
square optical lattice (lattice spacing ai a t = 532 nm) with 
variable lattice depths in x and y directions [23] (see Ap- 
pendix). A microscope objective with a resolution com- 
parable to the lattice spacing was used for fluorescence 
detection of individual atoms. Because inelastic light- 
assisted collisions during the imaging lead to a rapid loss 
of atom pairs, our scheme detects the parity of the atom 
number. We used an algorithm to deconvolve the im- 
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FIG. 1. Quantum-correlated particle-hole pairs 
in one-dimensional Mott insulators. (A) In a two- 
dimensional array of atoms (blue circles), decoupled one- 
dimensional systems were created by suppressing tunneling 
along the y direction. Quantum fluctuations then only induce 
correlated particle-hole excitations (yellow ellipses) along the 
x direction of the one-dimensional Mott insulators. (B) Such 
correlated fluctuations in the occupation hj are detected in 
the experiment as correlated fluctuations in the parity Sj . The 
light red bar in (A) marks the one-dimensional chain chosen 
in (B) for further explanations. 



ages, yielding single-site-resolved information of the on- 
site parity. Typically, our samples contained 150 — 200 
atoms in order to avoid Mis of occupation numbers n > 1. 

To detect particle-hole pairs, we evaluated two-site 
parity correlation functions |24j 



C(d) = (s k s k+d ) - (h){h+d), 



(2) 



where §k = e l7Tdnk is the parity operator at site k and 



d is the distance between the lattice sites. For the case 
of n = 1, Sfc yields +1(— 1) for an odd(even) occupation 
number n k . If a particle-hole pair exists on sites k and 
k + d, the same parity s(rik) = s(rik+d) = — 1 is detected 
(Fig. [I]). The existence of correlated particle- hole pairs 
therefore leads to an increase of (sfcSfc+d) above the fac- 
torized form (sk)(sk+d), which results from uncorrelated 
fluctuations, e.g., due to thermal excitations. We ob- 
tained C(d) from our deconvolved images by an average 
over many experimental realizations and by an additional 
average over k in a central region of interest. 

We first analyzed two-site parity correlations in one- 
dimensional systems (Fig.[2j\ and B). To create iso- 
lated one-dimensional tubes, we kept the lattice axis 
along y at a constant depth of V y = 17(1) EV, where 



E r = h 2 / (8maj 2 at ) denotes the recoil energy and m is the 
atomic mass of 87 Rb. We recorded the nearest-neighbor 
correlations C(d = 1) for different values of J/U along 
the direction of the one-dimensional tubes (red circles in 
Fig. [2j3) , where J and U are the tunneling matrix element 
and the on-site interaction energy in the Bosc-Hubbard 
model, respectively (see Appendix). For small J/U, the 
nearest-neighbor correlations vanish, since only uncorre- 
lated thermal excitations exist deep in the MI regime. 
As particle-hole pairs emerge with increasing J/U, we 
observe an increase of nearest-neighbor correlations un- 
til a peak value is reached, well before the critical value 
( J/U)l d « 0.3 IIU H6] for the one-dimensional SF-MI 
transition. The observed signal is a genuine quantum 
effect because thermally induced particle-hole pairs ex- 
tend over arbitrary distances and are therefore uncorre- 
lated. Their presence leads to a reduction of the corre- 
lation signal. We found no correlations when performing 
the same analysis perpendicular to the one-dimensional 
tubes (blue circles in Fig.[2j3) , showing that the coupling 
between the tubes was negligible. 

Our data show very good agreement with ab-initio 
finite-temperature Matrix Product State (MPS) calcu- 
lations gang at temperature T = 0.09 U /fee (Fig.[2^, 
solid line) that also take into account our harmonic trap- 
ping potential with frequency lu/(2it) = 60(l)Hz. Com- 
pared to a homogeneous system at T — (dashed- line) , 
the experimental signal is reduced, especially around the 
maximum. This reduction can be attributed in equal 
parts to the finite temperature of our system and the av- 
eraging over different local chemical potentials. The lat- 
ter is especially severe in the one-dimensional case owing 
to the narrow width of the Mott lobe for n = 1 close 
to the critical point [TS]. Interestingly, the growth of 
particle-hole correlations oc J 2 /U 2 expected from first- 
order perturbation theory (see Appendix) is limited to 
very small values J/U < 0.05, before significant devia- 
tions in the experiment and the numerical simulations 
are observed. 

As the dimensionality of the system plays an im- 
portant role in its correlation properties, we also mea- 
sured the two-site parity correlations across the two- 
dimensional SF-MI transition by simultaneously varying 
J/U along both lattice axes (Fig.[2]C). In contrast to the 
one-dimensional case, we now clearly observe the same 
nearest-neighbor correlations within our error bars along 
both axes. The maximum correlations are significantly 
smaller than in one dimension, and the peak value is 
now reached around the critical value [J/U) 2d ps 0.06 
[29] . We compared our data with Quantum Monte 
Carlo (QMC) simulations for a homogeneous system at 
T = 0.lU/k B (solid line in Fig.[2p) and found good 
quantitative agreement. Here, the broader shape of the 
Mott lobe leads to a weaker averaging effect over dif- 
ferent local chemical potentials. The increased strength 
of the correlations and the larger shift of the maximum 
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FIG. 2. Non-local parity correlations. (A) Top row: 
Typical experimental fluorescence images for J/U = 0.06 
(Al), J/U = 0.11 (A2) and J/U = 0.3 (A3) for the one- 
dimensional geometry. Bottom row: Reconstructed on-site 
parity. Particle-hole pairs are emphasized by a yellow shad- 
ing in (Al). For increased J/U the pairs start to proliferate 
and an identification in a single experimental image becomes 
impossible (A2, A3) (B) One-dimensional nearest-neighbor 
correlations C(d = 1) as a function of J/U along the x 
(red circles) and y directions (blue circles). The curves are 
first-order perturbation theory (dashed-dotted line), Density- 
Matrix Renormalization Group (DMRG) calculations for a 
homogeneous system at T — (dashed line) and finite- 
temperature MPS calculations including harmonic confine- 
ment at T = 0.09U/kB (solid line). (C) Parity correlations 
in two dimensions. Symbols have the same meaning as in (B). 
The curves are first-order perturbation theory (dashed-dotted 
line) and a QMC calculation for a homogeneous system at 
T = 0.01U/k B (dashed line) and T = Q.lU/k B (solid line). 
Each data point is an average over the central 9x7 lattice 
sites from 50 — 100 pictures. The error bars denote the la 
statistical uncertainty. The light blue shading highlights the 
SF phase. 



of the correlations relative to the critical point in the 
one-dimensional case directly reflect the more prominent 
role of quantum fluctuations in lower dimensions. This 
can also be seen from the on-site fluctuations C(d = 0) 
at the critical point which are significantly increased in 
the one-dimensional case (see Appendix). In both the 
one-dimensional and the two-dimensional systems, two- 
site correlations are expected to decay strongly with dis- 
tance. Our data for the next-nearest-neighbor correla- 
tions C(d = 2) is consistent with this predicted behavior 
(see Appendix). 

In addition to two-site correlations, we evaluated 
string-type correlators O p{ l) = (jljtl *,) (see Eq.gand 
Fig.jlb) in our one-dimensional systems, where the prod- 
uct is calculated over a chain of length I + 1. In the 
simplest case of a zero-temperature MI at J/U — 0, no 
fluctuations exist and therefore Op (I) = 1. As J/U in- 
creases, fluctuations in the form of particle-hole pairs ap- 
pear. Whenever a certain number of particle-hole pairs 
lies completely within the region covered by the string 
correlator, the respective minus signs cancel pairwisc. 
However, there is also the possibility that a particle-hole 
pair is cut by one end of the string correlator, for ex- 
ample when a particle exists at position < k and the 
corresponding hole has a position > k, resulting in an 
unpaired minus sign. As a consequence, Op(l) decreases 
with increasing J/U since the probability to cut particle- 
hole pairs becomes larger. Finally, at the transition to 
the SF phase, the pairs begin to deconfine and overlap, 
resulting in a completely random product of signs and 
Op = 0. 

To support this intuitive argument, we calculated 
Op(l) numerically using DMRG in a homogeneous sys- 
tem at T = 0. We show 0%(l) in Fig.[3]for selected dis- 
tances I together with the extrapolated values of Op — 
limz-j.oo Op(l) (inset to Fig 3), which we computed using 
finite-size scaling (see Appendix). We performed a fit to 
the extrapolated values close to the critical point with an 
exponential scaling 



2 pKe^(~A[{J/U)l d -{J/U)] 1/2 ) 



(3) 



characteristic for a transition of BKT type (Fig.[3| that 
one expects in one dimension [141 115] . From the fit we 
find (J/U) 1 / = 0.295 - 0.320, which is compatible with 
previously computed values [55J US] (see Appendix) . The 
fact that Op = in the SF and O p > in the MI as well 
as the agreement with the expected scaling show that Op 
serves as an order parameter for the MI phase in one di- 
mension Additionally, the simulations demonstrate 
that Op (I) is well suited to characterize the SF-MI tran- 
sition even for finite lengths I. 

Our experimentally obtained values of Op (I) for string 
length I < 8 (Fig.[4]A.) agree qualitatively well with in- 
trap MPS calculations at T = 0.09U/k B (Fig.g^). We 
observe a stronger decay of Op(l) with I compared to 
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FIG. 3. Numerical calculation of the string-order 
parameter. 0%(l) as a function of J/U calculated with 
DMRG for a homogeneous chain (n = 1, T = 0) of to- 
tal length 216. Lines show Op(l) for selected lengths I 
(black to red colors). Inset: Extrapolated value Op = 
lim^oo Op(l) together with a fit (black line) of the form 



1/2. 



Op oc exp ( - A [{J/U)l d - (J/U) j 
transition of BKT type (see Appendix 



characteristic for a 



the T = case, because at finite temperature thermal 
fluctuations lead to minus signs at random positions of 
the chain and reduce the average value of Op (I). Despite 
that, we still see a strong growth of Op (I), once the tran- 
sition from the SF to the MI is crossed, with a similar 
behavior as in Fig. [3] 

For a completely uncorrelated state, Op (I) factorizes 
t° Yik<j<k+i (*j)> an d m a homogeneous system we would 
expect a decay with string length of the form (§j) +1 , 
which can be slow provided the mean on-site parity (§j) 
is close to one. To rule out that our experimental data 
shows only such a trivial behavior, we define a new quan- 
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FIG. 4. String correlators. (A) Experimental values of 
Op (I) for lengths < I < 8 (B) In-trap MPS calculations 
at T = 0.09U/kp- (C) Experimentally determined string 
correlator Op(l) as defined in Eq.|i]for lengths 1 < I < 8 (D) 
In-trap MPS calculations at T = 0.09 U/k B . The I axes in 
(C) and (D) have been inverted. 



tity Op(l) that more naturally reflects the underlying 
correlations: 



Op 



(0 



op 



(0 



n 



(4) 



k<j<k+l 



First, we notice that Op(l) for length I = 1 is equal 
to the two-site correlation function C(d = 1). Second, 
Op (I) « Op (I) for long distances I since rife<j<fe+;(*i) 
eventually decays to zero (except for the singular case 
J/U = and T = 0). The correlation function Op (I) can 
therefore be understood as an extension of the two-site 
correlation function, that essentially captures the physics 
behind string order in one-dimensional Mis. 

Experimental and theoretical values for Op (I) are 
shown in Figs, and D. For small J/U, Op(l) is signifi- 
cantly reduced compared to Op(l) since few particle-hole 
pairs exist and Op(l) is close to its factorized form for 
short lengths I. In the case of vanishing J/U, we even 
expect Op(l) — since all sites are completely decou- 
pled. For intermediate J/U ~ 0.1, Op (I) grows rapidly 
with length / showing a strong deviation from the fac- 
torized form. Finally, in the SF regime, Op{l) becomes 
indiscernible from zero for large lengths in contrast to the 
nearest-neighbor two-site correlation function. Further- 
more, our data shows a genuine three-site correlation, 
which we revealed after subtracting all two site correla- 
tors in addition to local terms (see Appendix). 

We have shown direct measurements of non-local 
parity-parity correlation functions on the single-lattice- 
site and single-atom level and we demonstrated that a 
one-dimensional Mott insulator is characterized by non- 
local string order. A natural extension of our work would 
be to reveal, e.g., topological quantum phases such as 
the Haldane insulator of bosonic atoms jl3l [14] . A Hal- 
dane insulator exhibits a hidden antiferromagnetic order- 
ing and is expected to occur in one-dimensional quantum 
gases in the presence of longer ranged interactions, which 
could be realized in our experiment using Rydberg atoms 
(301. 
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APPENDIX 

Preparation of a two-dimensional degenerate 
quantum gas 

Most of the experimental details are the same as in 
Refs. [23] and [31], but several improvements have been 
implemented in order to increase the stability of the sys- 
tem. We started by transporting a thermal cloud of 
87 Rb atoms in the hyperfine state \F — 1, mp — —1) with 
a single-beam optical dipole trap (wavelength 1064 nm, 
beam waist Wg — 40 /im) in front of the high-resolution 
imaging system. Afterwards, the atoms were loaded 
into a crossed dipole trap, which consists of the hori- 
zontal lattice beams (wavelength 1064 nm, beam waist 
wq = 70 /im) whose retro-reflections were blocked with 
mechanical shutters. We evaporatively cooled the atoms 
in this trap by changing its depth from 20 /iK to 10 /iK 
within 1 s, before transferring the atoms into a ver- 
tical standing wave (wavelength 1064 nm, beam waist 
w = 70 ^m) of which we typically populate 40 antin- 
odes. 

To extract a single two-dimensional system, we used 
a position-dependent microwave transfer in a magnetic 
field gradient of dB/dz = 45 G/cm, 1.9 times larger com- 
pared to Ref. [23- After a transfer of all atoms from 
\F = 1, mjr = — 1) to \F = 2, nip = —2) with a broad mi- 
crowave sweep, the atoms from one antinode of the stand- 
ing wave were transferred back to \F — l,mp = —1) us- 
ing a HSl-pulse |31j of 5 ms duration and a sweep width 
of 3.5 kHz. The remaining atoms in state F = 2 were 
removed from the trap by a laser pulse resonant with the 
F = 2 ->• F' = 3 transition. 

For the final evaporation, we added a focussed laser 
beam of wavelength 850 nm shone in through the imaging 
system, with a beam diameter of ~ 20 fiin at the position 
of the atoms. The intensity of the beam is adjusted in 
such a way that a small number of atoms is lost due to 
three-body recombination, resulting in a stabilization of 
the atom number. In this configuration, we evaporatively 
cooled the atoms by tilting the trap horizontally with a 



G 



magnetic field gradient of 11.25 G/cm and by reducing 
the power of the 850 nm beam. We adjusted the atom 
number by changing the end point of this evaporation 
and achieved an atom number stability better than 20 %. 



Experimental sequence and lattice calibration 

For the measurements on the one-dimensional systems, 
we changed the potential depth of lattice axis y to a 
fixed final value V y = 17(1) E r using an s-shaped ramp of 
120 ms duration after having created a degenerate quan- 
tum gas in a single antinode of the vertical lattice beam. 
At the same time and using the same ramp shape, we 
ramped the lattice depth along x to variable final depths 
V x to attain different values of J/U. The vertical lattice 
depth was kept at V z — 20(1) E r during these ramps. For 
the measurement of two-site correlations in two dimen- 
sions, the depths of lattice axis x and y were changed 
simultaneously to the same final value. 

We obtained J and U from the lattice depths by a nu- 
merical band-structure calculation. In order to calibrate 
the lattice depths, we modulated the lattice intensity and 
measured parametric excitations of the atoms from the 
zcroth to the second Bloch band. As a result, we ob- 
served a resonance in the in-situ width of the atom cloud 
as a function of the modulation frequency. Our method 
allowed us to measure the transition frequency with an 
uncertainty of 2%. Besides this, we saw drifts of the 
lattice depth during our measurements of typically 5%. 
Both errors together lead to an experimental uncertainty 
of - 10% for J/U. 

For the detection, we froze the density distribution 
of the atoms by increasing all axes simultaneously to 
V x = 72 E r , V y = 78 E r , and V z = 83 E r within 0.2ms. 
We checked that our ramps were fast enough to avoid any 
local number-squeezing dynamics [22]. Finally, the indi- 
vidual atoms on the lattice were detected using a high- 
resolution objective with numerical aperture NA = 0.68. 
For this, we further increased the lattice depths along 
all three directions to ~ 3000 E r within 2 ms and illumi- 
nated the atoms with an optical molasses that induces 
fluorescence and simultaneously laser cools the atoms 



Bose-Hubbard model and perturbation theory 

For the temperatures and lattice parameters relevant 
to the presented measurements, the physics of a gas of 
interacting bosons in a lattice is captured by the Bose- 
Hubbard model [5]: 

H = -J ^ a\aj + y nt(ftj - 1) + ej-hj (5) 

<i,j> i i 



Here, a, and a\ correspond to the bosonic annihilation 
and creation operators at site i, hi = a\di is the on-site 
number operator, and Ci denotes the energy offset due to 
an external harmonic confinement. 

The emergence of quantum-correlated particle-hole 
pairs in a Mott insulator with finite tunnelling can be 
readily understood within first-order perturbation the- 
ory. For a Mott insulator in the atomic limit, where the 
ratio of the tunneling energy J to the on-site interac- 
tion energy U vanishes, the many-body state is simply 
given by \^f)j/u=o = FJi \ n i — 1) f° r the case °f unity 
filling and T = 0. For finite but still small tunnelling, 
J/U <C 1, the ground state can be approximated by con- 
sidering the tunnelling term as a perturbation. To first 
order, one obtains: 

\*)j/u<i oc \*)j/u=o + 7jY, «i & i\*)j/u=o, (6) 

(hi) 

which yields a probability to find a particle-hole pair 
on neighboring sites proportional to (J/U) 2 . Within 
the approximation of Eq.[6] the nearest-neighbor parity- 
correlation function is given by C(d = 1) = 16(J/C/) 2 + 
0[( J/U) 4 ]. Closer to the transition to the superfluid and 
for large system sizes, higher-order perturbation terms 
become more important and essentially lead to a rapid 
increase of bound particle-hole pairs and an extension of 
their size, eventually resulting in deconfinement of the 
pairs at the transition point. 

On-site variance and next-nearest-neighbor 
correlations 

In addition to the nearest-neighbor parity correlation 
C(d = 1), we also evaluated the correlations C(d) (Eq.[2| 
for distances d — and d = 2 in our one-dimensional 
systems (Fig.[5]Al). For d = 0, this amounts to measur- 
ing the on-site variance C(d = 0) = er(sfc) 2 . Deep in the 
Mott-insulating regime, the on-site number distribution 
is strongly squeezed and the variance is close to zero, 
whereas in the SF regime, it saturates at cr(sfc) 2 = 1, 
as expected for the case of a parity detection^] C(d) 
drops rapidly as a function of the distance d and the 
numerical calculations predict only a small maximum of 
0.01 in the next-nearest-neighbor correlation C(d = 2) 
at J/U ~ 0.17, which is however indiscernible from the 
statistical noise in our measurements (Fig.[5]A2). 

We also show the parity correlation C(d) in the two- 
dimensional system for different distances (Fig.[5|Bl). 



1 In Ref. 1231 we used a different definition of the parity operator, 
Pk = (s*; + l)/2, which yielded 1(0) for an odd (even) occupation 
number. In that case, the maximum variance was cr(p^) 2 = 
CT (s fc ) 2 /4 = 0.25. 
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FIG. 5. On-site variance and next-nearest-neighbor 
correlations. (Al) Parity correlations C(d) for different dis- 
tances d — 0,1,2 (green, red and gray circles, respectively) 
for the one-dimensional systems. The solid lines are finite- 
temperature MPS calculations including harmonic confine- 
ment at T — 0.09U/k B . (A2) Enlarged view of C(d) for 
d — 1,2 for the same datasets as (Al). (Bl) Same quantities 
as in (Al), but for the two-dimensional case whereas the data 
is the mean of C'(d) for both possible directions. (B2) En- 
larged view of C(d) for d = 1, 2 (red, gray) and additionally 
the two-site correlations for the next-nearest-neighbor along 
the diagonal (blue). The solid lines are a QMC calculation 
for a homogeneous system at T — 0.1(//fcs. The statisti- 
cal errorbars in (Al) and (Bl) are smaller than the dot size. 
We attribute the systematic shift of C(d = 0) in (B) to the 
inhomogeneous trapping potential. 



We found that, similar to the one-dimensional systems, 
next-nearest-neighbor correlations are not visible above 
the statistical noise, whereas the QMC simulations show 
a small maximum around J/U ~ 0.065 (Fig.|}B2). How- 
ever, a striking difference is that the on-site fluctuations 
around the critical point are significantly smaller than in 
one dimension. 



String order and multi-site correlations 

In the following paragraph, we show that our data for 
the string-type correlators (Fig. [4]) cannot be explained 
with pure two-site correlations. Additionally, we address 
the connection between string order and multi-site cor- 
relations. 

To illustrate this, let us consider the simplest case of 
a string- type correlator including three sites (S1S2S3), 
where si,§2 and S3 refer to the parity of three neigh- 
boring sites on a one-dimensional chain. If one of the 
sites is not correlated with the others, we can calculate 
(S1S2S3) from two-site and on-site terms. For instance, 
if site 3 is not correlated with sites 1 and 2, we have 
(•S1S2S3) = (S1S2) (,§3). This fact can be generally ex- 
pressed using a three-site cumulant (siS2S3)c, defined as 



-Ci t2 (s 3 ) - C 2 , 3 (si) - Ci )3 (s 2 ) 



(7.) 



with two-site correlation functions C+ j — (si§j) 



The cumulant is a measure of the correla- 



tions between all three sites as it vanishes if one of 
the sites is not correlated with the others. Addition- 
ally, if the cumulant vanishes, Eq.[7] can be written as 
(sis 2 s 3 ) = (si)(s 2 )(s 3 ) + Ci : 2(s 3 ) + C 2 , 3 (si) + Ci :3 (s 2 ) 
and we therefore have a situation where (S1S2S3) does 
not contain more information than two-site and on-site 
terms. 

Our experimental values for the three-site cumulant 
(siS2S3)c (Fig(6j) show a significant signal for J/U < 0.20, 
in quantitative agreement with a MPS calculation at 
T = 0.09U/ks including the harmonic confinement. 
Note that this is a substantial effect, as the peak value 
constitutes about 40% of the peak value of <D 2 P {1 = 2), 
where only on-site terms have been subtracted. We draw 
two conclusions from this: First, our data shows a three- 
site correlation beyond simple two-site correlations. Sec- 
ond, our data for string-type expectation values cannot 
be expressed in terms of two-site correlations alone. Let 
us explain the latter statement in more detail. It turns 
out that if a three-site expectation value cannot be ex- 
pressed via two-site terms, then expectation values in- 
cluding more than three sites cannot be expressed in this 
way either. The four-site term, e.g., can be written as: 



{S1S2S3S4) 



(S1S2S3S4) 

+ (SiS 2 S4) 



+ \SlS2S 3 ) c {S 4 ) + ( 
(S3) + (siS 3 S4)c(s 2 ) 



S2S3S4) 
+ C (2) 



(Si) 

(8) 



where contains only terms including two-site and on- 
site terms. Even if there were no correlations between 
four sites (siS2S3S4)c = 0, we still had to include the 
non-vanishing three-site cumulants. Therefore (S1S2S3S4) 
cannot be expressed via two-site and on-site terms. The 
argument can easily be extended to longer strings and 




J/U 



FIG. 6. Three-site correlator. Our experimental values 
for three site cumulant (siS2S3) c (green circles) show the ex- 
istence of three-site correlations in our system. The curves 
are DMRG calculations for a homogeneous system at T = 
(dashed line) and finite-temperature MPS calculations includ- 
ing harmonic confinement at T = 0.09(7/fcs (solid line). 
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we conclude that our signal for string-type correlators 
including more than two-sites cannot be explained with 
pure two-site correlations. 

In general, every string-type correlator can be ex- 
pressed as a sum of products of cumulants similar to 
Eq.|H] Since cumulants are a measure of multi-site cor- 
relations, such an expansion gives us information about 
the contribution of multi-site correlations to string order. 
In the atomic limit where J/U ~ 0, even two-site corre- 
lations are absent, but string order is present. In this 
case, the string signal is trivially dominated by on-site 
terms. For small but finite J/U values, three-site corre- 
lations almost vanish and the string signal is dominated 
by two-site correlations and on-site terms. In the range 
of J/U « 0.1 — 0.2, three-site correlations build up and 
also contribute to the string-order signal. Generally, it 
would be interesting to know how this relation extends 
when approaching the critical point. One particularly 
interesting question is whether close to the critical point 
multi-site correlations between an infinite number of sites 
exist in the thermodynamic limit. Such an analysis is dif- 
ficult theoretically and experimentally, which can already 
be seen from the expression for the four-site cumulant. 
This is however beyond the scope of this manuscript and 
a matter of further investigation. 

Instead, we continue with an interpretation of the 
present three-site correlations. One explanation for such 
a signal would be the existence of pairs extending over 
three sites. However, another possibility can explain 
such correlations: Assume we restrict ourselves to a sys- 
tem of only three sites, where the parity on each site 
can take the values Si = ±1, i = 1,2,3. The sys- 
tem is fully characterized by the probability p(si, s 2 , S3) 
of finding the parities Si,s 2 and S3 on sites 1,2 and 3. 
In a situation, where only nearest-neighbor pairs exist 
with probability p nn , we have £>(+,+,+) = 1 — 2p nra , 
p(-> +) = P(+, -) -) = Pnn and p{si, s 2 , s 3 ) = for 
all other cases. Interestingly, we find a non-vanishing 
three-site cumulant (s 1 s 2 s 3 ) c = 16p^ n — 32p^ n in this sit- 
uation. The three-site correlation remains present even 
if we consider the three sites as a subsystem of a longer 
chain. We conclude that a three-site correlation can arise 
from next-neighbor pairs alone, simply because site 1 is 
correlated with site 2 and site 2 is correlated with site 3. 
However, our signal extends far beyond the region where 
first order perturbation theory is valid and we can there- 
fore assume that particle-hole pairs with an extension of 
more than two sites as well as more complicated clusters 
play a significant role in our situation. 

Finally, we would like to point out an important dif- 
ference between two-site correlators and string correla- 
tors: Averaging over many experimental realizations of 
the system, a two-site correlator at distance d sums up 
correlations from pairs with exactly an extension of d but 
also anti-correlations from all pairs with a different size. 
In contrast, a string- type correlator of length I sums up 



positive contributions from all pairs which have a size 
less than I and lie within the string length. This feature 
makes the string-type correlator more suitable to study 
the Mott insulating phase. In particular, the phase tran- 
sition in one dimension is marked by the vanishing of the 
string correlators for long lengths. It is not possible to 
extract the same information using only two-site correla- 
tors. 



Comparison of string correlator with theory 

For a more detailed comparison of the experimentally 
obtained string correlations with theory (Fig.|4|, we show 
Op (I) for different lengths I = 1 (Fig.[7J red circles), I = 4 
(blue circles) and I — 8 (green circles) together with MPS 
calculations at T = 0.09 U/ks including the harmonic 
confinement. We observe a good qualitative agreement 
of the data with the numerical simulations. Systematic 
errors of the experimental data can arise from different 
mean atom numbers for different J/U, leading to differ- 
ent local chemical potentials. A systematic discrepancy 
between theory and experiment can result from a small 
mismatch of the trapping frequencies or the calibration 
of J and U. Additionally, the theoretical prediction is 
calculated at fixed temperature (in units of U), but the 
experimental data is approximately taken at a constant 
entropy. In the latter case, we expect the temperature 
to scale with U for small J/U values, but it is not clear 
that this scaling remains valid for higher J/U values. The 
temperature for the numerical simulation was chosen to 
yield the best agreement over the full range of J/U val- 
ues. 
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FIG. 7. String correlator - comparison with theory. 

String correlator Op (I) for different lengths I = 1 (red cir- 
cles), I — 4 (blue circles) and I — 8 (green circles). These 
curves are cuts of the three-dimensional representation of the 
data, as shown in Fig.[4] Solid lines correspond to finite- 
temperature MPS calculations including harmonic confine- 
ment at T = 0.09 U/kp with the same color coding as the 
experimental data. 
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DMRG and finite-size scaling 

We use a DMRG code with open boundary conditions, 
which implements the conservation of the number of par- 
ticles. The considered system sizes span between L = 168 
and L = 256. The number of retained states m has been 
chosen in order to have a truncation error smaller than 
10 _8 [33]. This results in a relative error for the shown 
data which is smaller than 10 -3 , as estimated via com- 
parison of data for different m between 185 and 235. 

The string correlator is strongly affected by the pres- 
ence of boundaries and finit- size effects. We introduce 
the notation P (J/U,l,L) to address the expectation 
value of the string correlator of length I obtained from 
the numerical simulation of a system of length L. The I 
sites are always taken in the center of the system to min- 
imize boundary effects. For the finite size scaling in the 
inset of Fig.[3j we analyzed correlators Op(J/U,aL,L) 
for a fixed fraction a of the total length L (I = aL). 
We extrapolated Op(J/U, a) = lim L ^oo Op(J/U, aL, L) 
using a scaling of the form Op(J/U,aL, L) = a + b/L r> 





a = 1/4 


a = 1/3 


a = 1/2 


Uh,(wh] = [0-23,0.37] 


0.303 


0.306 


0.296 


[(£)i,(£)a] = [0.24,0.36] 


0.310 


0.311 


0.299 


Uh,(£h] = [0.25,0.35] 


0.319 


0.318 


0.305 


[(jfh>(vh] = [0.26,0.34] 


0.321 


0.319 


0.311 


U)i,(£h] = [0.27,0.33] 


0.317 


0.314 


0.309 



TABLE I. Results for (J/U) c for different relative lengths a 
and fitting intervals [(^7)1,(^7)2]- The fitting error for a given 
a and (^) 2 ] is below 0.003. 



To determine (J/U) c , we fitted the extrapolated val- 
ues with Op cx exp ( - A[(J/U)l d - {J/U)Y 1/2 ). The 
results for (J/U) c appear to be strongly dependent on 
the fitting interval [(J/U)x, (J /U)2\ and a, as shown in 
Table [I] This large systematic error could be reduced us- 
ing more advanced finite size-scaling methods [351 137) . 
which is beyond the scope of this manuscript. In the 
inset of Fig.[3j we show the fit using a = 1/2 and 
[(J/[/)i,(J/C/) 2 ] = [0.25,0.35] . 

Finite-temperature MPS 

Thermal equilibrium states for a finite chain of q- 
dimensional systems can be approximated by MPS us- 



ing the techniques introduced in [37] HE]- For a chain 
of length L, the density matrix can be expressed as a 
vector in a Hilbert space of larger dimension, q 2L . We 
approximate this by a MPS of length 2L, where sites 
2k and 2k + 1 correspond to the k-th site in the physi- 
cal system. For inverse temperature f3 — l/(kpT), the 
thermal state (up to normalization) is formally identical 
to the imaginary time evolution of the identity operator 
p w e~ m = e-i H te~i H . To obtain the MPS repre- 
sentation of this operator, we use a second-order Suzuki- 
Trotter decomposition of the exponentials, and apply the 
imaginary time evolution to the initial state correspond- 
ing to the identity operator 27]. 

Using this method, we have simulated the thermal 
states of a Bose-Hubbard chain with length L = 30, 
for up to n = 4 particles per site. We used an aver- 
age over systems with different total atom numbers N, 
corresponding to the typical atom number distribution 
in the experiment, as well as an average over the central 
9 sites. The error of the numerical simulation originates 
from the Trotter step S and the bond dimension param- 
eter D, which controls the number of parameters in the 
MPS ansatz. Both errors can be reduced by running 
simulations with increasing (decreasing) D (5). We have 
used D = 20, S = 0.05, for which the estimated error is 
below 10~ 3 . 



Quantum Monte Carlo calculations 

Numerical results at finite temperature are obtained by 
using the Quantum Monte Carlo worm algorithm [38] in 
the implementation of Ref. [39] ■ This is an unbiased and 
statistically exact path integral Monte Carlo algorithm, 
formulated in continuous imaginary time, in which two 
worm operators (corresponding to a and operators) 
perform local updates in an extended configuration space 
and hereby directly sample the Green function. These 
updates lead to a fast decorrelation between the config- 
urations resulting in an integrated autocorrelation time 
of order unity. A bootstrap analysis shows that the rel- 
ative error is below 10 -3 for the presented data. It was 
previously demonstrated that experiments on the Bose- 
Hubbard model are in one-to-one agreement with such 
first principles simulations for realistic system sizes (on 
the order of a million atoms) and temperatures [40] . In 
our simulations, we assume that temperature scales as 
T ~ U which was shown to be a good assumption in this 
parameter regime |41j . 



